Aplicação de Survival Analysis no controle de Churn de assinaturas em Telecom

Introdução

Um dos assuntos mais recorrentes em qualquer tipo de serviço de assinatura é como reduzir o Churn, dado que conquistar novos clientes é bem mais difícil (e caro) do que manter os antigos.

Cerca de 70% das empresas sabem que é mais barato manter um cliente do que ter que ir atrás de um novo.

Fazendo uma analogia simples, os serviços de assinatura são como uma espécie de sangue na corrente sanguínea de uma empresa e uma interrupção de qualquer natureza prejudica todo o negócio, dado que esse é um modelo de negócio que se baseia na recorrência de tarifação e não no desenvolvimento, ou mesmo venda de outros produtos.

Modelos de negócios baseados no volume de pessoas que estão dispostas a terem uma cobrança recorrente o negócio fica bem mais complicado, dado que diferentemente de produtos que tem uma eslasticidade maior o fluxo de receita é extremamente sujeito aos sabores do mercado e dos clientes.

Dentro desse cenário, todas as empresas que tem o seu fluxo de receita baseado nesse tipo de business, saber quando um cliente entrará em uma situação de saída através do cancelamento do serviço (Churn) é fundamental para criar mecanismos de retenção mais efetivos, ou mesmo criação de réguas de contato com os clientes para evitar ou minimizar a chance de um cliente sair da base de dados.

Sendo assim, qualquer mecanismo ou mesmo esforço para minimizar esse efeito é de grande valia. Nos baseamos na teoria estatística buscar respostas para a pergunta: Como diminuir o Churn? Como identificar um potencial cliente que irá entrar em uma situação de Churn? Quais estratégias seguir para minimizar esse Churn? Quais réguas de comunicação com os clientes devemos ter para entender os motivos que estão fazendo um assinante cancelar o serviço e quais são as estratégias de customer winback possíveis nesse cenário?

E pra responder essa pergunta, fomos buscar as respostas na análise de sobrevivência dado que essa área da estatística é uma das que lidam melhor em termos de probabilidade de tempo de vida, seja de materiais (e.g. tempo de falha de algum sistema mecânico), no tempo de vida de pessoas propriamente ditas (e.g. dado uma determinada posologia qual é a estimativa de um paciente sobreviver a um câncer), e no nosso caso quanto tempo de vida um assinante tem até deixar cancelar a sua assinatura.

Análise de Sobrevivência

A análise de sobreviência é uma técnica estatístisca que foi desenvolvida na medicina e tem como principal finalidade estimar o tempo de sobrevivência ou tempo de morte de um determinado paciente dentro de um horizonte do tempo.

O estimador de Kaplan-Meier (1958) utiliza uma função de sobrevivência que leva em consideração uma divisão entre o número de observações que não falharam no tempo t pelo número total de observações no estudo em que cada intervalo de tempo tem-se o número de falhas/mortes/churn distintos bem como é calculado o risco de acordo com o número de individuos restantes no tempo subsequente.

Já o estimador Nelson-Aalen (1978) é um estimador que tem as mesmas características do Kaplan-Meier, com a dferença que esse estimador trabalha com uma função de sobrevivência que é a cumulative hazard rate function.

Os elementos fundamentais para caracterização de um estudo que envolve análise de sobrevivência são, o tempo inicial, escala de medida do intervalo de tempo e se o evento de churn ocorreu.

Os principais artigos são de Aalen (1978), Kaplan-Meier (1958) e Cox (1972).

Esse post não tem como principal objetivo dar algum tipo de introdução à survival analysis, dado que tem muitas referências na internet sobre o assunto e não há nada a ser acrescentado nesse sentido.

Como a análise de cohort, a análise de sobrevivência tem como principal característica ser um estudo de natureza logituginal, isto é, os seus resultados tem uma caracterisca de temporalidade seja em aspectos de retrospecção, quanto em termos de perspectivas, isso é, tem uma resposta tipicamente temporal para um detemrminado evento de interesse.

No caso o que vamos usar como forma de comparação amostral é o comportamento longituginal de acordo com determinadas caracteristicas de amostragens diferentes ao longo do tempo, e os fatores que influenciam no churn.

As diferenciações estão claras através das covariáveis, mas que devido a questões de NDA não vamos postar aqui, mas para fins de exemplificação e comparação vamos considerar que todas as covariaveis são semelhantes.

Podemos dizer que a análise de sobrevivência aplicada em um caso de telecom, pode ajudar ter uma estimativa em forma de probabilidade em relação ao tempo em que uma assinatura vai durar até o evento de churn (cancelamento) e dessa forma elaborar estratégias para evitar esse evento, dado que adquirir um novo cliente é mais caro do que manter um novo e entra totalmente dentro de uma estratégia de Customer Winback

No nosso caso o tempo de falha ou tempo de morte, como estamos falando de serviços de assinaturas, o nosso evento de interesse seria o churn, ou cancelamento da assinatura. Em outras palavras teríamos algo do tipo Time-to-Churn. Guardem esse termo.

Metodologia

Usamos dados de dois produtos antigos em que os dados foram anonimizados e aplicados um hash de embaralhamento uniforme (que obedece uma distribuição específica) nos atributos por questões de privacidade, que são:

  • id = Identificador do registro
  • product = produto
  • channel = canal no qual o cliente entrou na base de dados
  • free_user = flag que indica se o cliente entrou na base em gratuídade ou não
  • user_plan = se o usuário é pré-pago ou pós-pago
  • t = tempo que o assinante está na base de dados
  • c = informa se o evento de interesse (no caso o churn (cancelamento da assinatura) ocorreu ou não

Eliminamos o efeito de censura à esquerda eliminando casos de reativações, dado que queriamos entender a jornada do assinante como um todo sem nenhum tipo de viés relativo a questões de customer winback. Em relação à censura à direita temos que já se passaram alguns meses desde que essa base de dados foi extraída.

Um aspecto técnico importante a ser considerado é que esses dois produtos estão em categorias de comparabilidade, dado que sem isso nenhum tipo de caractericação seria nula.

No fim dessa implementação teremos uma tabela de vida de cada um desses produtos.

Implementação

Primeiramente vamos importar as bibilotecas: Pandas (para manipulação de dados), matplotlib (para a geração de gráficos), e lifelines para aplicação da análise de sobrevivência


In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import lifelines

Após realizar a importação das bibliotecas, vamos ajustar o tamanho das imagens para uma melhor visualização.


In [2]:
%pylab inline
pylab.rcParams['figure.figsize'] = (14, 9)


Populating the interactive namespace from numpy and matplotlib

Vamos realizar o upload da nossa base de dados criando um objecto chamado df e usando a classe read_csv do pandas.


In [3]:
df = pd.read_csv('https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/07%20-%20Survival/survival_data.csv')

Vamos checar a nossa base de dados


In [4]:
df.head()


Out[4]:
id product channel free_user user_plan t c
0 3315 B HH 1 0 22 0
1 2372 A FF 1 1 16 0
2 1098 B HH 1 1 22 0
3 2758 B HH 1 1 4 1
4 2377 A FF 1 1 29 0

Então como podemos ver temos as 7 variáveis na nossa base de dados.

Na sequência vamos importar a biblioteca do lifeniles, em especial o estimador de KaplanMaier


In [5]:
from lifelines import KaplanMeierFitter

kmf = KaplanMeierFitter()

Após realizar a importação da classe relativa ao estimador de Kaplan Meier no objeto kmf, vamos atribuir as nossas variáveis de tempo (T) e evento de interesse (C)


In [6]:
T = df["t"]

C = df["c"]

O que foi feito anteriormente é que buscamos no dataframe df o array t e atribuímos no objeto T, e buscamos o array da coluna c no dataframe e atribuímos no objeto C.

Agora vamos chamar o método fit usando esses dois objetos no snippet abaixo.


In [7]:
kmf.fit(T, event_observed=C )


Out[7]:
<lifelines.KaplanMeierFitter: fitted with 10000 observations, 6000 censored>

Objeto ajustado, vamos agora ver o gráfico relativo a esse objeto usando o estimador de Kaplan Meier.


In [8]:
kmf.survival_function_.plot()
plt.title('Survival function of Service Valued Add Products');
plt.ylabel('Probability of Living (%)')
plt.xlabel('Lifespan of the subscription (in days)')


Out[8]:
<matplotlib.text.Text at 0x101b24a90>

Como podemos ver no gráfico, temos algumas observações pertinentes, quando tratamos a probabilidade de sobreviência desses dois produtos no agregado que são:

  • Logo no primeiro dia há uma redução substâncial do tempo de sobrevivência da assinatura em aproximadamente 22%;
  • Há um decaímento quase que linear depois do quinto dia de assinatura;
  • Depois do dia número 30, a probabilidade de sobrevivência de uma assinatura é de aproximadamente de 50%. Em outras palavras: depois de 30 dias, metade dos novos assinantes já estarão fora da base de dados.

No entando, vamos plotar a mesma função de sobreviência considerando os intervalos de confiança estatística.


In [9]:
kmf.plot()
plt.title('Survival function of Service Valued Add Products - Confidence Interval in 85-95%');
plt.ylabel('Probability of Living (%)')
plt.xlabel('Lifespan of the subscription')


Out[9]:
<matplotlib.text.Text at 0x10ad8e0f0>

Contudo nesse modelo inicial temos duas limitações claras que são:

  • Os dados no agregado não dizem muito em relação à dinâmicas que podem estar na especificidade de alguns atributos/dimensões
  • Não são exploradas as dimensões (ou quebras) de acordo com os atributos que vieram na base de dados
  • Não há a divisão por produto

Para isso, vamos começar a entrar no detalhe em relação a cada uma das dimensões e ver o que cada uma tem de influência em relação à função de sobrevivência.

Vamos começar realizando a quebra pela dimensão que determina se o cliente entrou via gratruídade ou não (free_user).


In [10]:
ax = plt.subplot(111)

free = (df["free_user"] == 1)
kmf.fit(T[free], event_observed=C[free], label="Free Users")
kmf.plot(ax=ax, ci_force_lines=True)
kmf.fit(T[~free], event_observed=C[~free], label="Non-Free Users")
kmf.plot(ax=ax, ci_force_lines=True)
plt.ylim(0,1);
plt.title("Lifespans of different subscription types");
plt.ylabel('Probability of Living (%)')
plt.xlabel('Lifespan')


Out[10]:
<matplotlib.text.Text at 0x10ad8e908>

Este gráfico apresenta algumas informações importantes para os primeiros insights em relação a cada uma das curvas de sobreviência em relação ao tipo de gratuídade oferecida como fator de influência para o churn que são:

  • Os assinantes que entram como não gratuítos (i.e. não tem nenhum tipo de gratuidade inicial) após o 15o dia apresenta um decaímento brutal de mais de 40% da chance de sobreviência (tratando-se do intervalo de confiança);
  • Após o 15o dia os assinantes que não desfrutam de gratuidade tem a sua curva de sobrevivência em uma relativa estabilidade em torno de 60% na probabilidade de sobrevivência até o período censurado;
  • Ainda nos usuários sem gratuidade, dado o grau de variabilidade do intervalo de confiança podemos tirar como conclusão que muitos cancelamentos estão ocorrendo de forma muito acelerada, o que deve ser investigado com mais calma pelo time de produtos;
  • Já os usuários que entram via gratuidade (i.e. ganham alguns dias grátis antes de serem tarifados) apresenta um nível de decaímento do nível de sobreviência maior seja no período incial, quando ao longo do tempo, contudo uma estabilidade é encontrada ao longo de toda a série sem maiores sobressaltos.

Dado essa análise incial das curvas de sobreviência, vamos avaliar agora as probabilidades de sobrevivência de acordo com o produto.


In [11]:
ax = plt.subplot(111)

product = (df["product"] == "A")
kmf.fit(T[product], event_observed=C[product], label="Product A")
kmf.plot(ax=ax, ci_force_lines=True)
kmf.fit(T[~product], event_observed=C[~product], label="Product B")
kmf.plot(ax=ax, ci_force_lines=True)

plt.ylim(0,1);
plt.title("Survival Curves of different Products");
plt.ylabel('Probability of Living (%)')
plt.xlabel('Lifespan')


Out[11]:
<matplotlib.text.Text at 0x10aeaabe0>

Este gráfico apresenta a primeira distinção entre os dois produtos de uma forma mais clara.

Mesmo com os intervalos de confiança com uma variação de 5%, podemos ver que o produto A (linha azul) tem uma maior probabilidade de sobrevivência com uma diferença percentual de mais de 15%; diferença essa amplificada depois do vigésimo dia.

Em outras palavras: Dado um determinada safra de usuários, caso o usuário entre no produto A o mesmo tem uma probabilidade de retenção de cerca de 15% em relação a um usuário que por ventura entre no produto B, ou o produto A apresenta uma cauda de retenção superior ao produto B.

Empiricamente é sabido que um dos principais fatores de influência de produtos SVA são os canais de mídia os quais esses produtos são oferecidos.

O canal de mídia é o termômetro em que podemos saber se estamos oferencendo os nossos produtos para o público alvo correto.

No entanto para um melhor entendimento, vamos analisar os canais nos quais as assinaturas são originadas.

A priori vamos normalizar a variável channel para realizar a segmentação dos canais de acordo com o conjunto de dados.


In [12]:
df['channel'] = df['channel'].astype('category');
channels = df['channel'].unique()

Após normalização e transformação da variável para o tipo categórico, vamos ver como está o array.


In [13]:
channels


Out[13]:
[HH, FF, CC, AA, GG, ..., BB, EE, DD, JJ, ZZ]
Length: 11
Categories (11, object): [HH, FF, CC, AA, ..., EE, DD, JJ, ZZ]

Aqui temos a representação de 11 canais de mídia os quais os clientes entraram no serviço.

Com esses canais, vamos identificar a probabilidade de sobrevivência de acordo com o canal.


In [14]:
for i,channel_type in enumerate(channels):
    ax = plt.subplot(3,4,i+1)
    ix = df['channel'] == channel_type
    kmf.fit( T[ix], C[ix], label=channel_type )
    kmf.plot(ax=ax, legend=True)
    plt.title(channel_type)
    plt.xlim(0,40)
    if i==0:
        plt.ylabel('Probability of Survival by Channel (%)')
plt.tight_layout()


Fazendo uma análise sobre cada um desses gráficos temos algumas considerações sobre cada um dos canais:

  • HH, DD: Uma alta taxa de mortalidade (churn) logo antes dos primeiros 5 dias, o que indica uma caracteristica de efemeridade ou atratividade no produto para o público desse canal de mídia.
  • FF: Apresenta menos de 10% de taxa de mortalidade nos primeiros 20 dias, e tem um padrão muito particular depois do 25o dia em que praticamente não tem uma mortalidade tão alta. Contém um intervalo de confiança com uma oscilação muito forte.
  • CC: Junto com o HH apesar de ter uma taxa de mortalidade alta antes do 10o dia, apresenta um grau de previsibildiade muito bom, o que pode ser utilizado em estratégias de incentivos de mídia que tenham que ter uma segurança maior em termos de retenção a médio prazo.
  • GG, BB: Apresentam uma boa taxa de sobrevivência no inicio do período, contudo possuem oscilações severas em seus respectivos intervalos de confiança. Essa variável deve ser considerada no momento de elaboração de uma estratégia de investimento nesses canais.
  • JJ: Se houvesse uma definição de incerteza em termos de sobreviência, esse canal seria o seu melhor representante. Com os seus intervalos de confiança oscilando em mais de 40% em relação ao limite inferior e superior, esse canal de mídia mostra-se extremamente arriscado para os investimentos, dado que não há nenhum tipo de regularidade/previsibilidade de acordo com esses dados.
  • II: Apesar de ter um bom grau de previsibilidade em relação à taxa de sobreviência nos primeiros 10 dias, após esse período tem uma curva de hazard muito severa, o que indica que esse tipo de canal pode ser usado em uma estratégia de curto prazo.
  • AA, EE, ZZ: Por haver alguma forma de censura nos dados, necessitam de mais análise nesse primeiro momento. (Entrar no detalhe dos dados e ver se é censura à direita ou algum tipo de truncamento).

Agora que já sabemos um pouco da dinâmica de cada canal, vamos criar uma tabela de vida para esses dados.

A tabela de vida nada mais é do que uma representação da função de sobrevivencia de forma tabular em relação aos dias de sobrevivência. (Melhorar essa parte).

Para isso vamos usar a biblioteca utils do lifelines para chegarmos nesse valor.


In [15]:
from lifelines.utils import survival_table_from_events

Biblioteca importada, vamos usar agora as nossas váriaveis T e C novamente para realizar o ajuste da tabela de vida.


In [16]:
lifetable = survival_table_from_events(T, C)

Tabela importada, vamos dar uma olhada no conjunto de dados.


In [17]:
print (lifetable)


          removed  observed  censored  entrance  at_risk
event_at                                                
0            2250      2247         3     10000    10000
1             676       531       145         0     7750
2             482       337       145         0     7074
3             185       129        56         0     6592
4             232        94       138         0     6407
5             299        85       214         0     6175
6             191        73       118         0     5876
7             127        76        51         0     5685
8             211        75       136         0     5558
9            2924        21      2903         0     5347
10            121        27        94         0     2423
11             46        27        19         0     2302
12             78        26        52         0     2256
13            111        16        95         0     2178
14             55        35        20         0     2067
15            107        29        78         0     2012
16            286        30       256         0     1905
17            156        23       133         0     1619
18            108        18        90         0     1463
19             49        11        38         0     1355
20             50        17        33         0     1306
21             61        13        48         0     1256
22            236        23       213         0     1195
23             99         6        93         0      959
24            168         9       159         0      860
25            171         7       164         0      692
26             58         6        52         0      521
27             77         2        75         0      463
28             29         6        23         0      386
29            105         1       104         0      357
30             69         0        69         0      252
31            183         0       183         0      183

Diferentemente do R que possuí a tabela de vida com a porcentagem relativa à probabilidade de sobrevivência, nesse caso vamos ter que fazer um pequeno ajuste para obter a porcentagem de acordo com o atributo entrance e at_risk.

O ajuste se dará da seguinte forma:


In [18]:
survivaltable = lifetable.at_risk/np.amax(lifetable.entrance)

Ajustes efetuados, vamos ver como está a nossa tabela de vida.


In [19]:
survivaltable


Out[19]:
event_at
0     1.0000
1     0.7750
2     0.7074
3     0.6592
4     0.6407
5     0.6175
6     0.5876
7     0.5685
8     0.5558
9     0.5347
10    0.2423
11    0.2302
12    0.2256
13    0.2178
14    0.2067
15    0.2012
16    0.1905
17    0.1619
18    0.1463
19    0.1355
20    0.1306
21    0.1256
22    0.1195
23    0.0959
24    0.0860
25    0.0692
26    0.0521
27    0.0463
28    0.0386
29    0.0357
30    0.0252
31    0.0183
Name: at_risk, dtype: float64

Vamos transformar a nossa tabela de vida em um objeto do pandas para melhor manipulação do conjunto de dados.


In [20]:
survtable = pd.DataFrame(survivaltable)

Para casos de atualização de Churn-at-Risk podemos definir uma função que já terá a tabela de vida e poderá fazer a atribuição da probabilidade de sobreviência de acordo com os dias de sobreviência.

Para isso vamos fazer uma função simples usando o próprio python.


In [21]:
def survival_probability( int ):
   survtable["at_risk"].iloc[int]
   print ("The probability of Survival after", int, "days is", survtable["at_risk"].iloc[int]*100, "%") 
   return;

Nesse caso vamos ver a chance de sobreviência usando o nosso modelo Kaplan-Meier já ajustado para uma assinatura que tenha 22 dias de vida.


In [22]:
survival_probability(22)


The probability of Survival after 22 days is 11.95 %

Ou seja, essa assinatura tem apenas 11.95% de probabilidade de estar ativa, o que significa que em algum momento muito próximo ela pode vir a ser cancelada.

Conclusão

Como podemos ver acima, usando análise de sobrevivência podemos tirar insights interessantes em relação ao nosso conjunto de dados, em especial para descobrirmos a duração das assinaturas em nossa base de dados, e estimar um tempo até o evento de churn.

Os dados utilizados refletem o comportamento de dois produtos reais, porém, que foram anonimizados por questões óbvias de NDA. Contudo nada impede a utilização e a adaptação desse código para outros experimentos. Um ponto importante em relação a essa base de dados é que como pode ser observado temos uma censura à direita muito acentuada o que limita um pouco a visão dos dados a longo prazo, principalmente se houver algum tipo de cauda longa no evento de churn.

Como coloquei no [São Paulo Big Data Meetup de Março] (http://www.slideshare.net/esportesocial/sp-big-data-meetup-march16) há uma série de arquiteturas que podem ser combinadas com esse tipo de análise, em especial métodos de Deep Learning que podem ser um endpoint de um pipeline de predição.

Espero que tenham gostado e quaisquer dúvidas mandem uma mensagem para flavioclesio at gmail.com

PS: Agradecimentos especiais aos meus colegas e revisores Eiti Kimura, Gabriel Franco e Fernanda Eleuterio.